library(ggplot2)
library(scales)

# 1.加载所需库------
if (!require("pacman")) install.packages("pacman")
pkgs <- c("dplyr", "ggplot2", "readxl","writexl", "tidyr", "lubridate", "janitor", "stringr", "purrr",
          "psych", "magrittr","scales") #, "sf", "spdep", "tmap", "ggspatial")
pacman::p_load(pkgs, character.only = TRUE)


source("./code/2.3 人口数估发病率.R")

# 读取区县代码
geo_code <- read_excel("data/origin/st/GEO_副本.xlsx") %>% 
  dplyr::select(LocationIDs = `location ID`,county_code = new_location_id) %>% 
  distinct()

st <- readxl::read_excel("./data/origin/st_result/时空扫描-2005-2023-逐年.xlsx", 
              sheet = "Sheet2") %>% filter(`P-value` <= 0.05) 

names(st)
st_clean <- st %>%
  dplyr::select(year, LocationIDs = `Location IDs`,Cluster_type) %>% 
  mutate(LocationIDs = str_replace_all(LocationIDs, "\\s+", "")) %>%
  separate_rows(LocationIDs, sep = ",") %>%
  filter(!is.na(LocationIDs)) %>%
  arrange(year, LocationIDs) %>% left_join(geo_code) %>% ungroup() %>% 
  dplyr::select(-LocationIDs)

writexl::write_xlsx(st_clean,"./data/origin/st_result/st_clean.xlsx")

st_count <- st_clean %>% group_by(year,Cluster_type) %>% 
  summarise(count = n()) 


# 数据预处理
plot_data <- st_count %>%
  mutate(Cluster_type = factor(Cluster_type,
                               levels = c("cluster II", "cluster I"),
                               labels = c("Secondary Cluster", "Primary Cluster")))

# 期刊推荐配色 (ColorBrewer Paired)
cluster_colors <- c("Primary Cluster" = "#3182BD", 
                    "Secondary Cluster" = "#6BAED6")

# 创建图形
cluster_plot <- ggplot(plot_data, 
                       aes(x = factor(year), 
                           y = count, 
                           fill = Cluster_type)) +
  geom_col(position = "stack", 
           width = 0.7, 
           color = "white", 
           linewidth = 0.3) +
  scale_fill_manual(values = cluster_colors) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + # 优化Y轴间距
  labs(title = "Spatiotemporal Clustering of Brucellosis in Xinjiang (2005-2023)",
       x = "Year",
       y = "Number of Clustered Counties",
       fill = "Cluster Type") +
  theme_minimal(base_size = 12) +
  theme(
    text = element_text(family = "Arial"),
    plot.title = element_text(size = 14, hjust = 0.5, face = "bold",
                              margin = margin(b = 15)),
    axis.title = element_text(size = 12),
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1,
                               color = "gray30"),
    axis.text.y = element_text(color = "gray30"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "top",
    legend.justification = "right",
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 11),
    plot.margin = unit(c(1,1,1,1), "cm")
  )

cluster_plot

# 高质量输出
ggsave("./主图/时空扫描地图/聚类柱状堆叠图.tiff", 
       plot = cluster_plot,
       device = "tiff",
       dpi = 600,
       width = 25,   # 适合双栏显示
       height = 15,  # 单位：cm
       units = "cm",
       compression = "lzw")


st_table <- st %>% 
  mutate(count = str_count(`Location IDs`, ",") + 1) %>% 
  dplyr::select(`Scan time frame` = year, `Cluster time` = `Time frame`, 
                `Centroid (latitude, longitude)/radius (km)` = `Coordinates / radius`, 
                `Cluster counties (n)` = count, `LLR` = "Log likelihood ratio", `RR` = "Relative risk", `P-value`,p_numeric = `P-value`) %>% 
  mutate(
    Significance = case_when(
      p_numeric < 0.001                ~ "极显著 (P < 0.001)",
      p_numeric >= 0.001 & p_numeric <= 0.05 ~ "显著 (0.001 ≤ P ≤ 0.05)",
      p_numeric > 0.05                 ~ "不显著 (P > 0.05)",
      TRUE                              ~ NA_character_
    ) %>% ordered(levels = c("极显著 (P < 0.001)", 
                             "显著 (0.001 ≤ P ≤ 0.05)",
                             "不显著 (P > 0.05)"))
  ) %>% 
  filter(Significance %in% c("极显著 (P < 0.001)", "显著 (0.001 ≤ P ≤ 0.05)"))




writexl::write_xlsx(st_table,"./data/origin/st_result/聚类逐年个数.xlsx")


# # 假设你的数据框名为st_clean
# st_clean1 <- st_clean %>% 
#   left_join(nameCE %>% filter(region == "county") %>% select(county_code, nameC)) %>% 
#   # 按照year和Cluster_type分组
#   group_by(year, Cluster_type) %>%
#   # 对每个分组，将nameC用逗号连接，并新增一列merged_name
#   summarise(county_name = paste(nameC, collapse = ",")) %>%
#   # 重新排列列的顺序
#   select(year, county_name, Cluster_type)


# # 数据清洗流程
# cleaned_st <- st %>%
#   # 预处理：移除Location IDs中的空格，确保分隔符统一
#   mutate(`Location IDs` = str_replace_all(`Location IDs`, "\\s+", "")) %>%
#   # 拆分多值Location IDs为独立行
#   separate_rows(`Location IDs`, sep = ",") %>%
#   # 转换Location IDs为数值类型
#   # mutate(`Location IDs` = as.numeric(`Location IDs`)) %>%
#   # 移除可能产生的空值
#   filter(!is.na(`Location IDs`)) %>%
#   # 按原始列顺序排序
#   select(names(st)[1], `Location IDs`, everything())
